#Intake is a package to load and share datasets, and will be our connection to the cloud via the
#intake-esm catalogue.
import intake
#Xarray is essentially Pandas for n-dimensional datasets (which will be the outputs from
#climate models)
import xarray as xr
#Proplot is the 'next big thing' for data visualisation in Python, apparently.
import proplot
#NetCDF for importing observed datasets
import netCDF4
#These modules are imported for obvious reasons
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set(rc={'figure.facecolor':'white'})
plt.style.use('default')
#Parallelisation
from joblib import Parallel, delayed
#remove any warning messages
import warnings
warnings.filterwarnings('ignore')
We read the data catalogue provided by https://pangeo.io/cloud.html:
#necessary url
url = "https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json"
#open the catalog
dataframe = intake.open_esm_datastore(url)
Check that this dataframe has actually worked:
dataframe
pangeo-cmip6 catalog with 7686 dataset(s) from 514961 asset(s):
| unique | |
|---|---|
| activity_id | 18 |
| institution_id | 36 |
| source_id | 88 |
| experiment_id | 170 |
| member_id | 657 |
| table_id | 37 |
| variable_id | 700 |
| grid_label | 10 |
| zstore | 514961 |
| dcpp_init_year | 60 |
| version | 737 |
We can see this tells us the number of different model catagories we are able to choose from.
We now choose the model from this list, with the following specification:
#Historical and AMIP model ensembles
model_ensembles = dataframe.search(source_id='GISS-E2-1-G',
experiment_id=['historical', 'amip'],
table_id='Amon',
variable_id=['tas','pr','psl'])
This query yields an intake_esm.core.esm_datastore data type, which we can use in order to get the dataset we have searched for.
Before we do this, let's learn a bit more about the model ensembles we're importing:
model_ensembles
pangeo-cmip6 catalog with 2 dataset(s) from 198 asset(s):
| unique | |
|---|---|
| activity_id | 1 |
| institution_id | 1 |
| source_id | 1 |
| experiment_id | 2 |
| member_id | 47 |
| table_id | 1 |
| variable_id | 3 |
| grid_label | 1 |
| zstore | 198 |
| dcpp_init_year | 0 |
| version | 12 |
print(model_ensembles.df.head().to_latex(index=False))
\begin{tabular}{lllllllllrr}
\toprule
activity\_id & institution\_id & source\_id & experiment\_id & member\_id & table\_id & variable\_id & grid\_label & zstore & dcpp\_init\_year & version \\
\midrule
CMIP & NASA-GISS & GISS-E2-1-G & historical & r2i1p1f1 & Amon & tas & gn & gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/his... & NaN & 20180827 \\
CMIP & NASA-GISS & GISS-E2-1-G & historical & r2i1p1f1 & Amon & pr & gn & gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/his... & NaN & 20180827 \\
CMIP & NASA-GISS & GISS-E2-1-G & historical & r1i1p1f1 & Amon & tas & gn & gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/his... & NaN & 20180827 \\
CMIP & NASA-GISS & GISS-E2-1-G & historical & r2i1p1f1 & Amon & psl & gn & gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/his... & NaN & 20180827 \\
CMIP & NASA-GISS & GISS-E2-1-G & historical & r1i1p1f1 & Amon & psl & gn & gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/his... & NaN & 20180827 \\
\bottomrule
\end{tabular}
model_ensembles.df.loc[model_ensembles.df['experiment_id'] == 'historical'].head()
| activity_id | institution_id | source_id | experiment_id | member_id | table_id | variable_id | grid_label | zstore | dcpp_init_year | version | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | CMIP | NASA-GISS | GISS-E2-1-G | historical | r2i1p1f1 | Amon | tas | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/hi... | NaN | 20180827 |
| 1 | CMIP | NASA-GISS | GISS-E2-1-G | historical | r2i1p1f1 | Amon | pr | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/hi... | NaN | 20180827 |
| 2 | CMIP | NASA-GISS | GISS-E2-1-G | historical | r1i1p1f1 | Amon | tas | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/hi... | NaN | 20180827 |
| 3 | CMIP | NASA-GISS | GISS-E2-1-G | historical | r2i1p1f1 | Amon | psl | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/hi... | NaN | 20180827 |
| 4 | CMIP | NASA-GISS | GISS-E2-1-G | historical | r1i1p1f1 | Amon | psl | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/hi... | NaN | 20180827 |
print(
'The number of simulations in the historical ensemble for temperature at surface (tas) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'historical')
& (model_ensembles.df['variable_id'] == 'tas')].shape[0])
The number of simulations in the historical ensemble for temperature at surface (tas) is: 46
print(
'The number of simulations in the historical ensemble for precipitation flux (pr)) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'historical')
& (model_ensembles.df['variable_id'] == 'pr')].shape[0])
The number of simulations in the historical ensemble for precipitation flux (pr)) is: 46
print(
'The number of simulations in the historical ensemble for air pressure at sea level (psl) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'historical')
& (model_ensembles.df['variable_id'] == 'psl')].shape[0])
The number of simulations in the historical ensemble for air pressure at sea level (psl) is: 46
Notice: The number of historical simulations for each variable is the same, which is good news!
model_ensembles.df.loc[model_ensembles.df['experiment_id'] == 'amip'].head()
| activity_id | institution_id | source_id | experiment_id | member_id | table_id | variable_id | grid_label | zstore | dcpp_init_year | version | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 30 | CMIP | NASA-GISS | GISS-E2-1-G | amip | r5i1p1f1 | Amon | pr | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/am... | NaN | 20181016 |
| 31 | CMIP | NASA-GISS | GISS-E2-1-G | amip | r5i1p1f1 | Amon | tas | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/am... | NaN | 20181016 |
| 32 | CMIP | NASA-GISS | GISS-E2-1-G | amip | r4i1p1f1 | Amon | psl | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/am... | NaN | 20181016 |
| 33 | CMIP | NASA-GISS | GISS-E2-1-G | amip | r5i1p1f1 | Amon | psl | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/am... | NaN | 20181016 |
| 34 | CMIP | NASA-GISS | GISS-E2-1-G | amip | r4i1p1f1 | Amon | tas | gn | gs://cmip6/CMIP6/CMIP/NASA-GISS/GISS-E2-1-G/am... | NaN | 20181016 |
print(
'The number of simulations in the AMIP ensemble for temperature at surface (tas) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'amip')
& (model_ensembles.df['variable_id'] == 'tas')].shape[0])
The number of simulations in the AMIP ensemble for temperature at surface (tas) is: 20
print(
'The number of simulations in the AMIP ensemble for precipitation flux (pr) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'amip')
& (model_ensembles.df['variable_id'] == 'pr')].shape[0])
The number of simulations in the AMIP ensemble for precipitation flux (pr) is: 20
print(
'The number of simulations in the AMIP ensemble for air pressure at sea level (psl) is:',
model_ensembles.df.loc[
(model_ensembles.df['experiment_id'] == 'amip')
& (model_ensembles.df['variable_id'] == 'psl')].shape[0])
The number of simulations in the AMIP ensemble for air pressure at sea level (psl) is: 20
Notice: The number of AMIP simulations for each variable is the same, which is good news!
model_ensembles_dataset = model_ensembles.to_dataset_dict()
--> The keys in the returned dictionary of datasets are constructed as follows: 'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'
We check what we have downloaded:
model_ensembles_dataset.keys()
dict_keys(['CMIP.NASA-GISS.GISS-E2-1-G.amip.Amon.gn', 'CMIP.NASA-GISS.GISS-E2-1-G.historical.Amon.gn'])
#AMIP datasets
model_amip_ds = model_ensembles_dataset['CMIP.NASA-GISS.GISS-E2-1-G.amip.Amon.gn']
#Historical datasets
model_historical_ds = model_ensembles_dataset['CMIP.NASA-GISS.GISS-E2-1-G.historical.Amon.gn']
The above are xarray.core.dataset.Dataset straight away, meaning they can be readily used for anything one would want with the Xarray package (which is suited to work with gridded meteorological data).
model_amip_ds
<xarray.Dataset>
Dimensions: (member_id: 20, lat: 90, bnds: 2, lon: 144, time: 1980)
Coordinates:
* member_id (member_id) <U8 'r1i1p1f1' 'r1i1p1f2' ... 'r5i1p3f1' 'r5i1p5f1'
* lat (lat) float64 -89.0 -87.0 -85.0 -83.0 ... 83.0 85.0 87.0 89.0
lat_bnds (lat, bnds) float64 dask.array<chunksize=(90, 2), meta=np.ndarray>
* lon (lon) float64 1.25 3.75 6.25 8.75 ... 351.2 353.8 356.2 358.8
lon_bnds (lon, bnds) float64 dask.array<chunksize=(144, 2), meta=np.ndarray>
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
time_bnds (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
height float64 ...
Dimensions without coordinates: bnds
Data variables:
pr (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 600, 90, 144), meta=np.ndarray>
psl (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 963, 90, 144), meta=np.ndarray>
tas (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 963, 90, 144), meta=np.ndarray>
Attributes: (12/40)
title: GISS-E2-1-G output prepared for CMIP6
contact: Kenneth Lo (cdkkl@giss.nasa.gov)
institution_id: NASA-GISS
variable_id: tas
grid_label: gn
branch_time_in_child: 0.0
... ...
source: GISS-E2.1G (2019): \naerosol: Varies with physic...
grid: atmospheric grid: 144x90
external_variables: areacella
source_id: GISS-E2-1-G
sub_experiment_id: none
intake_esm_dataset_key: CMIP.NASA-GISS.GISS-E2-1-G.amip.Amon.gnWe verify that this has lat x lon = 90 x 144, 20 members in the ensemble, and spans 1850-2014.
Moreover, the variables tas, pr and psl are float32 values within an array (ensemble, time, lat, lon).
model_historical_ds
<xarray.Dataset>
Dimensions: (member_id: 46, lat: 90, bnds: 2, lon: 144, time: 1980)
Coordinates:
* member_id (member_id) <U10 'r101i1p1f1' 'r102i1p1f1' ... 'r9i1p5f1'
* lat (lat) float64 -89.0 -87.0 -85.0 -83.0 ... 83.0 85.0 87.0 89.0
lat_bnds (lat, bnds) float64 dask.array<chunksize=(90, 2), meta=np.ndarray>
* lon (lon) float64 1.25 3.75 6.25 8.75 ... 351.2 353.8 356.2 358.8
lon_bnds (lon, bnds) float64 dask.array<chunksize=(144, 2), meta=np.ndarray>
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
time_bnds (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
height float64 ...
Dimensions without coordinates: bnds
Data variables:
pr (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 600, 90, 144), meta=np.ndarray>
psl (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 600, 90, 144), meta=np.ndarray>
tas (member_id, time, lat, lon) float32 dask.array<chunksize=(1, 600, 90, 144), meta=np.ndarray>
Attributes: (12/46)
title: GISS-E2-1-G output prepared for CMIP6
contact: Kenneth Lo (cdkkl@giss.nasa.gov)
institution_id: NASA-GISS
parent_experiment_id: piControl
variable_id: tas
grid_label: gn
... ...
source: GISS-E2.1G (2019): \naerosol: Varies with physic...
grid: atmospheric grid: 144x90, ocean grid: 288x180
external_variables: areacella
source_id: GISS-E2-1-G
sub_experiment_id: none
intake_esm_dataset_key: CMIP.NASA-GISS.GISS-E2-1-G.historical.Amon.gnWe verify that this has lat x lon = 90 x 144, 46 members in the ensemble, and spans 1850-2014.
Moreover, the variables tas, pr and psl are float32 values within an array (ensemble, time, lat, lon).
This gives a nice geometric understanding of the tas array:
model_historical_ds.tas
<xarray.DataArray 'tas' (member_id: 46, time: 1980, lat: 90, lon: 144)>
dask.array<getitem, shape=(46, 1980, 90, 144), dtype=float32, chunksize=(1, 600, 90, 144), chunktype=numpy.ndarray>
Coordinates:
* member_id (member_id) <U10 'r101i1p1f1' 'r102i1p1f1' ... 'r9i1p5f1'
* lat (lat) float64 -89.0 -87.0 -85.0 -83.0 ... 83.0 85.0 87.0 89.0
* lon (lon) float64 1.25 3.75 6.25 8.75 ... 351.2 353.8 356.2 358.8
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
height float64 ...
Attributes:
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
long_name: Near-Surface Air Temperature
standard_name: air_temperature
units: KNot needed for this model.
Due to different griddings of each model, we have to 'guess' where Australia should be and verify on a projection plot:
#Set lat and lon indices for each experiment
A_lat = range(22, 41)
A_lon = range(45,65)
#Tas variable
model_historical_A_tas_ds = model_historical_ds['tas'][:, :, A_lat, A_lon]
model_amip_A_tas_ds = model_amip_ds['tas'][:, :, A_lat, A_lon]
#Pr variable
model_historical_A_pr_ds = model_historical_ds['pr'][:, :, A_lat, A_lon]
model_amip_A_pr_ds = model_amip_ds['pr'][:, :, A_lat, A_lon]
#Psl variable
model_historical_A_psl_ds = model_historical_ds['psl'][:, :, A_lat, A_lon]
model_amip_A_psl_ds = model_amip_ds['psl'][:, :, A_lat, A_lon]
model_historical_A_tas_ds
<xarray.DataArray 'tas' (member_id: 46, time: 1980, lat: 19, lon: 20)>
dask.array<getitem, shape=(46, 1980, 19, 20), dtype=float32, chunksize=(1, 600, 19, 20), chunktype=numpy.ndarray>
Coordinates:
* member_id (member_id) <U10 'r101i1p1f1' 'r102i1p1f1' ... 'r9i1p5f1'
* lat (lat) float64 -45.0 -43.0 -41.0 -39.0 ... -15.0 -13.0 -11.0 -9.0
* lon (lon) float64 113.8 116.2 118.8 121.2 ... 153.8 156.2 158.8 161.2
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
height float64 ...
Attributes:
cell_measures: area: areacella
cell_methods: area: time: mean
comment: near-surface (usually, 2 meter) air temperature
long_name: Near-Surface Air Temperature
standard_name: air_temperature
units: Kmodel_historical_A_pr_ds
<xarray.DataArray 'pr' (member_id: 46, time: 1980, lat: 19, lon: 20)>
dask.array<getitem, shape=(46, 1980, 19, 20), dtype=float32, chunksize=(1, 600, 19, 20), chunktype=numpy.ndarray>
Coordinates:
* member_id (member_id) <U10 'r101i1p1f1' 'r102i1p1f1' ... 'r9i1p5f1'
* lat (lat) float64 -45.0 -43.0 -41.0 -39.0 ... -15.0 -13.0 -11.0 -9.0
* lon (lon) float64 113.8 116.2 118.8 121.2 ... 153.8 156.2 158.8 161.2
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
height float64 ...
Attributes:
cell_measures: area: areacella
cell_methods: area: time: mean
comment: includes both liquid and solid phases
long_name: Precipitation
standard_name: precipitation_flux
units: kg m-2 s-1model_historical_A_psl_ds
<xarray.DataArray 'psl' (member_id: 46, time: 1980, lat: 19, lon: 20)>
dask.array<getitem, shape=(46, 1980, 19, 20), dtype=float32, chunksize=(1, 600, 19, 20), chunktype=numpy.ndarray>
Coordinates:
* member_id (member_id) <U10 'r101i1p1f1' 'r102i1p1f1' ... 'r9i1p5f1'
* lat (lat) float64 -45.0 -43.0 -41.0 -39.0 ... -15.0 -13.0 -11.0 -9.0
* lon (lon) float64 113.8 116.2 118.8 121.2 ... 153.8 156.2 158.8 161.2
* time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
height float64 ...
Attributes:
cell_measures: area: areacella
cell_methods: area: time: mean
comment: Sea Level Pressure
long_name: Sea Level Pressure
standard_name: air_pressure_at_sea_level
units: PaWe can see that these datasets are restricted to the geographical location of Ireland.
We create a plot of the first few members of the historical ensemble in the first month:
We plot some ensembles:
#Global historical Near-Surface Temperature Ensembles in 1850-01
model_historical_ds.tas.isel(time=0, member_id=range(0, 6)).plot(col='member_id',
col_wrap=3,
robust=True)
plt.suptitle("GISS Global Distribution of Tas in 1850-01 in Historical Ensemble", y=1.1)
plt.savefig('images/Global_historical_tas_ensemble', bbox_inches='tight', dpi=300)
#Global historical Precipitation Flux Ensembles in 1850-01
model_historical_ds.pr.isel(time=0, member_id=range(0, 6)).plot(col='member_id',
col_wrap=3,
robust=True)
plt.suptitle("GISS Global Distribution of Pr in 1850-01 in Historical Ensemble", y=1.1)
plt.savefig('images/Global_historical_pr_ensemble', bbox_inches='tight', dpi=300)
#Global historical Air Pressure at Sea Level Ensembles in 1850-01
model_historical_ds.psl.isel(time=0, member_id=range(0, 6)).plot(col='member_id',
col_wrap=3,
robust=True)
plt.suptitle("GISS Global Distribution of Psl in 1850-01 in Historical Ensemble", y=1.1)
plt.savefig('images/Global_historical_psl_ensemble', bbox_inches='tight', dpi=300)
We now create some Robinson plots for the first ensemble member of each variable:
#Create a Robinson projection with the proplot package
fig, ax = proplot.subplots(
axwidth=4.5,
tight=True,
proj='robin',
proj_kw={'lon_0': 180},
)
#Format options
ax.format(
land=False,
coast=True,
innerborders=True,
borders=True,
labels=True,
geogridlinewidth=0,
)
#Add contour variables
map1 = ax.contourf(model_historical_ds['lon'],
model_historical_ds['lat'],
#Choose ensemble member 0 and time 0
model_historical_ds['tas'][0, 0, :, :],
cmap='viridis',
extend='both')
#Add a colorbar below
ax.colorbar(map1, loc='r', shrink=1, extendrect=True)
plt.show()
#Create a Robinson projection with the proplot package
fig, ax = proplot.subplots(
axwidth=4.5,
tight=True,
proj='robin',
proj_kw={'lon_0': 180},
)
#Format options
ax.format(
land=False,
coast=True,
innerborders=True,
borders=True,
labels=True,
geogridlinewidth=0,
)
#Add contour variables
map1 = ax.contourf(model_historical_ds['lon'],
model_historical_ds['lat'],
#Choose ensemble member 0 and time 0
model_historical_ds['pr'][0, 0, :, :],
cmap='viridis',
extend='both')
#Add a colorbar below
ax.colorbar(map1, loc='r', shrink=1, extendrect=True)
plt.show()
plt.savefig('images/Global_historical_pr_time0_member1', bbox_inches='tight', dpi=300)
<Figure size 640x480 with 0 Axes>
#Create a Robinson projection with the proplot package
fig, ax = proplot.subplots(
axwidth=4.5,
tight=True,
proj='robin',
proj_kw={'lon_0': 180},
)
#Format options
ax.format(
land=False,
coast=True,
innerborders=True,
borders=True,
labels=True,
geogridlinewidth=0,
)
#Add contour variables
map1 = ax.contourf(model_historical_ds['lon'],
model_historical_ds['lat'],
#Choose ensemble member 0 and time 0
model_historical_ds['psl'][0, 0, :, :],
cmap='viridis',
extend='both')
#Add a colorbar below
ax.colorbar(map1, loc='r', shrink=1, extendrect=True)
plt.show()
plt.savefig('images/Global_historical_psl_time0_member1', bbox_inches='tight', dpi=300)
<Figure size 640x480 with 0 Axes>
...and also an Australia based one:
#ensure the map is high quality
proplot.rc.reso = 'hi'
#create a Robinson projection with the proplot package
fig, ax = proplot.subplots(
axwidth=4.5,
tight=True,
proj='robin',
proj_kw={'lon_0': 0},
)
#format options
ax.format(
land=False,
coast=True,
innerborders=True,
borders=True,
labels=True,
geogridlinewidth=0,
)
ax[0].format(latlim=(-43, -7), lonlim=(110, 160), labels=True)
#add contour variables
map1 = ax.contourf(model_historical_A_tas_ds['lon'],
model_historical_A_tas_ds['lat'],
model_historical_A_tas_ds[0, 0, :, :],
cmap='IceFire',
extend='both')
#add a colorbar below
ax.colorbar(map1, loc='r', shrink=1, extendrect=True)
plt.show()
plt.savefig('images/Au_historical_tas', bbox_inches='tight', dpi=300)
<Figure size 640x480 with 0 Axes>
#Australia historical Near-Surface Temperature Ensembles in 1850-01
model_historical_A_tas_ds.isel(time=0, member_id=range(0, 6)).plot(col='member_id',
col_wrap=3,
robust=True)
plt.suptitle("GISS Australia Distribution of Tas in 1850-01 in Historical Ensemble", y=1.1)
plt.savefig('images/Au_historical_tas_ensemble', bbox_inches='tight', dpi=300)
We will use the following algorithm to generate a probability distribution for the tas variable in each experiment:
(1): Match the time-span that each experiment covers.
(2): Across the first ten years and the last ten years of the simulation, calculate the global mean tas for each member of the ensemble.
(3): Calculate the difference in this mean.
(4): Normalise the data and plot as a histogram.
We create a function that calculates the mean difference in tas, or any other variable for that matter, given the starting time, final time, number of members in the ensemble, and the experiment type.
As this code normally takes a long time to run, we use parallelisation to speed this up. We also save the output array as a text file that we can load in after:
def change_mean(model, initialTime, totalTime, totalMembers, experiment, variable):
#Calculate 10 years in months
firstTenYears = range(initialTime, initialTime+12*10)
lastTenYears = range(totalTime-12*10,totalTime)
def _calculate_means(member):
first_mean_var = np.mean(model[member, firstTenYears, :, :])
last_mean_var = np.mean(model[member, lastTenYears, :, :])
return first_mean_var, last_mean_var
#first_mean_var = model[member, firstTenYears, :, :].mean(dim=['time','lat','lon'])
#last_mean_var = model[member, lastTenYears, :, :].mean(dim=['time','lat','lon'])
# Get the results
results = Parallel(n_jobs=-1)(
delayed(_calculate_means)(member) for member in range(totalMembers)
)
results = np.array(results)
first_mean_var = results[:, 0]
last_mean_var = results[:, 1]
change_mean_var = last_mean_var - first_mean_var
#Save the array as a txt file:
np.savetxt('arrays/change_mean_'+experiment+'_'+variable+'.txt', change_mean_var)
#Return the array as well
return change_mean_var
We calculate this for the historical ensemble.
totalTime = len(model_historical_A_tas_ds[0,:,0,0])
totalHistMembers = len(model_historical_A_tas_ds[:,0,0,0])
%%time
change_mean_historical_tas = change_mean(model=model_historical_A_tas_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalHistMembers,
experiment='historical',
variable='tas')
Run the following code if you have already run the function above:
change_mean_historical_tas = np.loadtxt('arrays/change_mean_historical_tas.txt')
Now we have an array of the mean temperature change in each ensemble member:
change_mean_historical_tas
array([0.8828125 , 0.97427368, 0.94189453, 0.84942627, 0.72814941,
1.17230225, 1.08071899, 0.6388855 , 0.69671631, 0.90206909,
0.99581909, 1.06787109, 0.849823 , 0.84683228, 0.77612305,
1.12915039, 0.91577148, 0.92495728, 0.81787109, 0.8241272 ,
1.01092529, 1.00643921, 0.93328857, 0.93548584, 0.96490479,
0.90756226, 1.2880249 , 0.84805298, 0.99835205, 0.83291626,
1.15353394, 0.76242065, 0.9541626 , 1.00869751, 1.1461792 ,
0.84570312, 0.89877319, 0.9850769 , 1.10644531, 0.86801147,
0.96920776, 1.13931274, 0.86953735, 0.85388184, 1.07110596,
1.07785034])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_historical_tas, hist=False, kde=True, rug=True, color='darkblue',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Tas (K)')
plt.title('GISS Historical Ensemble Distribution')
plt.savefig('images/Historical_ensemble_tas_dist', bbox_inches='tight', dpi=300)
We repeat this process above for each of the other two variables pr and psl:
%%time
change_mean_historical_pr = change_mean(model=model_historical_A_pr_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalHistMembers,
experiment='historical',
variable='pr')
change_mean_historical_pr = np.loadtxt('arrays/change_mean_historical_pr.txt')
change_mean_historical_pr
array([-1.00323450e-06, -8.85960617e-07, 9.73925125e-07, -2.14140528e-07,
2.04390562e-06, 7.05058483e-08, -1.06584775e-06, -4.25794497e-07,
1.68062616e-06, -9.34287527e-07, 6.93629772e-07, -1.84176315e-07,
1.24959479e-06, -1.98332964e-06, -1.22145138e-07, 1.77800757e-07,
-8.61982699e-07, -1.56050919e-06, 2.31240847e-07, 9.07901267e-07,
4.85146302e-07, 4.60404408e-07, 1.18032403e-06, -1.26381019e-06,
-6.65493644e-07, 3.10730502e-06, 4.39704309e-07, 7.76220986e-07,
1.45675585e-07, 2.62947651e-07, -3.51355993e-08, 4.03451850e-08,
4.64378900e-07, 1.69555096e-06, 3.21715925e-06, 1.79009658e-06,
-1.19237302e-06, -6.45863111e-07, 1.12302041e-06, 1.04702303e-06,
2.94742676e-06, 2.11438601e-06, 2.16277840e-06, -2.32217462e-06,
1.21092125e-06, 2.89401214e-07])
We notice that these values are particularly small due to the units of $kg m^{-2}s^{-1}$, which will cause problems in the future. So, we scale up these values to $g m^{-2}d^{-1}$, which equates to scaling up by a value of: $$ 1000 * (24 * 60^2) = 8.64 * 10^{7} $$
change_mean_historical_pr *= (8.63e7)
change_mean_historical_pr
array([ -86.5791375 , -76.45840124, 84.04973833, -18.48032753,
176.38905538, 6.08465471, -91.98266125, -36.74606505,
145.03803723, -80.62901361, 59.86024935, -15.89441599,
107.84003025, -171.16134786, -10.54112545, 15.34420535,
-74.38910689, -134.67194312, 19.95608509, 78.35187935,
41.86812585, 39.7329004 , 101.86196414, -109.06681946,
-57.43210149, 268.16042318, 37.94648183, 66.98787111,
12.57180302, 22.6923823 , -3.03220222, 3.48178946,
40.07589905, 146.32604816, 277.64084298, 154.4853345 ,
-102.90179162, -55.73798644, 96.91666146, 90.3580878 ,
254.36292981, 182.47151238, 186.64777599, -200.40366999,
104.50250429, 24.97532478])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_historical_pr, hist=False, kde=True, rug=True, color='darkblue',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Pr ($g m^{-2} d^{-1}$)')
plt.title('GISS Historical Ensemble Distribution')
plt.savefig('images/Historical_ensemble_pr_dist', bbox_inches='tight', dpi=300)
%%time
change_mean_historical_psl = change_mean(model=model_historical_A_psl_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalHistMembers,
experiment='historical',
variable='psl')
change_mean_historical_psl = np.loadtxt('arrays/change_mean_historical_psl.txt')
change_mean_historical_psl
array([-32.09375 , -3.1875 , 4.4921875, -1.1484375, 27.140625 ,
24.0859375, 17.921875 , 22.0703125, -4.0703125, 3.921875 ,
21.0390625, -3.09375 , 20.4921875, 16.2578125, -0.796875 ,
8.828125 , 12.484375 , 26.90625 , -21.1484375, 20.078125 ,
-5.328125 , 38.6484375, 19.4921875, -6.6015625, 23.6015625,
-25.1484375, 10.21875 , -5.3359375, 25.3125 , 11.53125 ,
-23.046875 , -35.2265625, 2.75 , -19.6640625, 20.5078125,
33.9453125, 20.6953125, 28.6875 , 11.5546875, 26.34375 ,
-9.2109375, 20.9453125, 16.0546875, 7.453125 , -3.578125 ,
5.6796875])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_historical_psl, hist=False, kde=True, rug=True, color='darkblue',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Psl (Pa)')
plt.title('GISS Historical Ensemble Distribution')
plt.savefig('images/Historical_ensemble_psl_dist', bbox_inches='tight', dpi=300)
We repeat this for the AMIP ensemble.
totalTime = len(model_amip_A_tas_ds[0,:,0,0])
totalAMIPMembers = len(model_amip_A_tas_ds[:,0,0,0])
%%time
change_mean_amip_tas = change_mean(model=model_amip_A_tas_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalAMIPMembers,
experiment='amip',
variable='tas')
Run the following code if you have already run the function above:
change_mean_amip_tas = np.loadtxt('arrays/change_mean_amip_tas.txt')
Now we have an array of the mean temperature change in each ensemble member:
change_mean_amip_tas
array([0.9503479 , 0.947052 , 0.81710815, 0.77191162, 0.77172852,
0.84463501, 0.85787964, 0.82977295, 0.88287354, 0.81079102,
0.90890503, 0.82131958, 0.7885437 , 1.00582886, 0.96640015,
0.87503052, 0.89993286, 0.76248169, 0.78479004, 0.76168823])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_amip_tas, hist=False, kde=True, rug=True, color='darkred',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Tas (K)')
plt.title('GISS AMIP Ensemble Distribution')
plt.savefig('images/AMIP_ensemble_tas_dist', bbox_inches='tight', dpi=300)
%%time
change_mean_amip_pr = change_mean(model=model_amip_A_pr_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalAMIPMembers,
experiment='amip',
variable='pr')
change_mean_amip_pr = np.loadtxt('arrays/change_mean_amip_pr.txt')
change_mean_amip_pr
array([1.24086000e-06, 2.43743671e-06, 4.01226862e-06, 3.23656241e-06,
2.45798947e-06, 2.15924410e-06, 2.64887967e-06, 3.89296110e-06,
2.64395749e-06, 4.40221993e-06, 1.33689537e-06, 2.12137638e-06,
3.43019383e-06, 2.07390258e-06, 2.50316589e-06, 2.54348379e-06,
3.15585930e-06, 2.73275691e-06, 1.97992631e-06, 5.28448436e-06])
We notice that these values are particularly small due to the units of $kg m^{-2}s^{-1}$, which will cause problems in the future. So, we scale up these values to $g m^{-2}d^{-1}$, which equates to scaling up by a value of: $$ 1000 * (24 * 60^2) = 8.64 * 10^{7} $$
change_mean_amip_pr *= (8.63e7)
change_mean_amip_pr
array([107.08621812, 210.35078771, 346.25878179, 279.31533568,
212.12449101, 186.34276621, 228.5983157 , 335.96254325,
228.1735311 , 379.91157988, 115.37407008, 183.07478185,
296.0257274 , 178.97779253, 216.02321613, 219.50265091,
272.35065791, 235.83692146, 170.86764055, 456.05100022])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_amip_pr, hist=False, kde=True, rug=True, color='darkred',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Pr ($g m^{-2} d^{-1}$)')
plt.title('GISS AMIP Ensemble Distribution')
plt.savefig('images/AMIP_ensemble_pr_dist', bbox_inches='tight', dpi=300)
%%time
change_mean_amip_psl = change_mean(model=model_amip_A_psl_ds,
initialTime=0,
totalTime=totalTime,
totalMembers=totalAMIPMembers,
experiment='amip',
variable='psl')
change_mean_amip_psl = np.loadtxt('arrays/change_mean_amip_psl.txt')
change_mean_amip_psl
array([ -37.4375 , -46.453125 , -50.9609375, -10.8203125,
-38.2109375, 17.5390625, -13.203125 , -21.515625 ,
-64.125 , -52.6640625, -43.8046875, -10.4609375,
-70.21875 , -63.6015625, -117.90625 , -17.1328125,
-67.96875 , -23.7578125, -58.2578125, -34.484375 ])
#Create a rug plot
fig, ax = plt.subplots(figsize=(8, 6))
sns.distplot(change_mean_amip_psl, hist=False, kde=True, rug=True, color='darkred',
kde_kws={'linewidth': 3}, rug_kws={'color': 'black'})
plt.xlabel('Mean Change in Psl (Pa)')
plt.title('GISS AMIP Ensemble Distribution')
plt.savefig('images/AMIP_ensemble_psl_dist', bbox_inches='tight', dpi=300)
#Mathematical and statistical packages
import math
import scipy
#KL
from scipy.stats import entropy
#JS
from scipy.spatial.distance import jensenshannon
#Wasserstein
import ot
import ot.plot
#Energy
import dcor
We re-apply the process above, but for all three variables at once. Hence we need to create the two three-variate array distributions for a distance to be calculated:
historical_outputs = np.dstack([change_mean_historical_tas, change_mean_historical_pr, change_mean_historical_psl])[0,:,:]
historical_outputs.shape
(46, 3)
amip_outputs = np.dstack([change_mean_amip_tas, change_mean_amip_pr, change_mean_amip_psl])[0,:,:]
amip_outputs.shape
(20, 3)
We see that these both have shape (number of amip ensembles, number of variables).
np.savetxt('arrays/all_historical_outputs.txt', historical_outputs)
historical_outputs = np.loadtxt('arrays/all_historical_outputs.txt')
historical_outputs
array([[ 0.8828125 , -86.5791375 , -32.09375 ],
[ 0.97427368, -76.45840124, -3.1875 ],
[ 0.94189453, 84.04973833, 4.4921875 ],
[ 0.84942627, -18.48032753, -1.1484375 ],
[ 0.72814941, 176.38905538, 27.140625 ],
[ 1.17230225, 6.08465471, 24.0859375 ],
[ 1.08071899, -91.98266125, 17.921875 ],
[ 0.6388855 , -36.74606505, 22.0703125 ],
[ 0.69671631, 145.03803723, -4.0703125 ],
[ 0.90206909, -80.62901361, 3.921875 ],
[ 0.99581909, 59.86024935, 21.0390625 ],
[ 1.06787109, -15.89441599, -3.09375 ],
[ 0.849823 , 107.84003025, 20.4921875 ],
[ 0.84683228, -171.16134786, 16.2578125 ],
[ 0.77612305, -10.54112545, -0.796875 ],
[ 1.12915039, 15.34420535, 8.828125 ],
[ 0.91577148, -74.38910689, 12.484375 ],
[ 0.92495728, -134.67194312, 26.90625 ],
[ 0.81787109, 19.95608509, -21.1484375 ],
[ 0.8241272 , 78.35187935, 20.078125 ],
[ 1.01092529, 41.86812585, -5.328125 ],
[ 1.00643921, 39.7329004 , 38.6484375 ],
[ 0.93328857, 101.86196414, 19.4921875 ],
[ 0.93548584, -109.06681946, -6.6015625 ],
[ 0.96490479, -57.43210149, 23.6015625 ],
[ 0.90756226, 268.16042318, -25.1484375 ],
[ 1.2880249 , 37.94648183, 10.21875 ],
[ 0.84805298, 66.98787111, -5.3359375 ],
[ 0.99835205, 12.57180302, 25.3125 ],
[ 0.83291626, 22.6923823 , 11.53125 ],
[ 1.15353394, -3.03220222, -23.046875 ],
[ 0.76242065, 3.48178946, -35.2265625 ],
[ 0.9541626 , 40.07589905, 2.75 ],
[ 1.00869751, 146.32604816, -19.6640625 ],
[ 1.1461792 , 277.64084298, 20.5078125 ],
[ 0.84570312, 154.4853345 , 33.9453125 ],
[ 0.89877319, -102.90179162, 20.6953125 ],
[ 0.9850769 , -55.73798644, 28.6875 ],
[ 1.10644531, 96.91666146, 11.5546875 ],
[ 0.86801147, 90.3580878 , 26.34375 ],
[ 0.96920776, 254.36292981, -9.2109375 ],
[ 1.13931274, 182.47151238, 20.9453125 ],
[ 0.86953735, 186.64777599, 16.0546875 ],
[ 0.85388184, -200.40366999, 7.453125 ],
[ 1.07110596, 104.50250429, -3.578125 ],
[ 1.07785034, 24.97532478, 5.6796875 ]])
np.savetxt('arrays/all_amip_outputs.txt', amip_outputs)
amip_outputs = np.loadtxt('arrays/all_amip_outputs.txt')
amip_outputs
array([[ 0.9503479 , 107.08621812, -37.4375 ],
[ 0.947052 , 210.35078771, -46.453125 ],
[ 0.81710815, 346.25878179, -50.9609375 ],
[ 0.77191162, 279.31533568, -10.8203125 ],
[ 0.77172852, 212.12449101, -38.2109375 ],
[ 0.84463501, 186.34276621, 17.5390625 ],
[ 0.85787964, 228.5983157 , -13.203125 ],
[ 0.82977295, 335.96254325, -21.515625 ],
[ 0.88287354, 228.1735311 , -64.125 ],
[ 0.81079102, 379.91157988, -52.6640625 ],
[ 0.90890503, 115.37407008, -43.8046875 ],
[ 0.82131958, 183.07478185, -10.4609375 ],
[ 0.7885437 , 296.0257274 , -70.21875 ],
[ 1.00582886, 178.97779253, -63.6015625 ],
[ 0.96640015, 216.02321613, -117.90625 ],
[ 0.87503052, 219.50265091, -17.1328125 ],
[ 0.89993286, 272.35065791, -67.96875 ],
[ 0.76248169, 235.83692146, -23.7578125 ],
[ 0.78479004, 170.86764055, -58.2578125 ],
[ 0.76168823, 456.05100022, -34.484375 ]])
Let's see these outputs in a 3D projection plot:
fig = plt.figure(figsize=(10,8))
#Make 3D projection
ax = plt.axes(projection='3d')
ax.scatter3D(historical_outputs[:,0], historical_outputs[:,1], historical_outputs[:,2], s=50, cmap='Greens', edgecolor='black', label='Historical Outputs')
ax.scatter3D(amip_outputs[:,0], amip_outputs[:,1], amip_outputs[:,2], s=50, cmap='Reds', marker= '^', edgecolor='black', label='AMIP Outputs')
ax.set_xlabel('Tas (K)')
ax.set_ylabel('Pr ($g m^{-2} d^{-1}$)')
ax.set_zlabel('Psl (Pa)')
ax.set_title('GISS Distribution of Mean Change', fontsize='x-large')
plt.legend()
plt.savefig('images/Tas_pr_prt_dist.eps', format='eps', bbox_inches='tight', dpi=300)
The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.
This will be necessary in order to calculate KL divergence (and JS distance) in a sensible way, as unlike Wasserstein and energy distance, KL divergence requires a 1D pdf for each measure, with a shared support (this is one of the limitations of KL divergence!).
from scipy.stats import multivariate_normal
amip_mean = np.mean(amip_outputs, axis=0)
amip_cov = np.cov(amip_outputs, rowvar=0)
amip_mean
array([ 0.85295105, 242.91044047, -41.27226563])
amip_cov
array([[ 5.49843069e-03, -3.31286009e+00, -8.74872201e-01],
[-3.31286009e+00, 7.53017700e+03, -9.61951382e+01],
[-8.74872201e-01, -9.61951382e+01, 8.66467133e+02]])
historical_mean = np.mean(historical_outputs, axis=0)
historical_cov = np.cov(historical_outputs, rowvar=0)
historical_mean
array([ 0.94459666, 33.06244523, 7.61871603])
historical_cov
array([[ 1.83231723e-02, 1.34800646e+00, 1.62243091e-01],
[ 1.34800646e+00, 1.24771428e+04, -4.33603802e+01],
[ 1.62243091e-01, -4.33603802e+01, 3.16340960e+02]])
amip_gaussian = multivariate_normal.pdf(amip_outputs, mean=amip_mean, cov=amip_cov)
This gives us a 1D array of probabilities which can be used for KL divergence calculation:
amip_gaussian
array([1.01558483e-04, 1.69041306e-04, 2.12160122e-04, 1.99791991e-04,
1.05934256e-04, 4.53103344e-05, 2.48294826e-04, 1.51176047e-04,
3.25623737e-04, 1.21142294e-04, 1.50457261e-04, 1.88034768e-04,
1.03186498e-04, 4.88981384e-05, 1.30603364e-05, 2.49280498e-04,
2.40335466e-04, 1.49771981e-04, 2.64746817e-05, 2.03393766e-05])
We repeat this process for the historical outputs:
historical_gaussian = multivariate_normal.pdf(historical_outputs,
mean=historical_mean, cov=historical_cov)
historical_gaussian
array([1.02754569e-05, 1.13318450e-04, 2.10977754e-04, 1.54277572e-04,
1.08966888e-05, 3.87010319e-05, 6.11782602e-05, 1.07379701e-05,
1.97490877e-05, 1.34916327e-04, 1.65084653e-04, 1.08742437e-04,
1.03640563e-04, 3.39108743e-05, 9.72970444e-05, 9.01445844e-05,
1.43410907e-04, 4.41070954e-05, 4.54315148e-05, 1.05131906e-04,
1.57578118e-04, 4.88612299e-05, 1.53202754e-04, 7.47092678e-05,
1.14591812e-04, 4.87632403e-06, 9.34901548e-06, 1.37404011e-04,
1.34666833e-04, 1.62771930e-04, 1.18082769e-05, 6.32004571e-06,
2.27752292e-04, 4.01596342e-05, 7.35987310e-06, 2.73556642e-05,
8.47258563e-05, 8.40234555e-05, 1.03050051e-04, 9.32409753e-05,
2.21567917e-05, 3.18954380e-05, 6.29046567e-05, 2.37707478e-05,
1.04366540e-04, 1.42192766e-04])
As the probability arrays have different lengths, we randomly select the 'matching' number of elements from the historical probabilities array, calculate the KL divergence, then repeat many times in order to calculate the average KL score.
def uneven_arrays(p,q):
#Match the array size
if len(p) > len(q):
p = np.random.choice(p, len(q))
elif len(q) > len(p):
q = np.random.choice(q, len(p))
#Return the outputs
return p, q
def uneven_KL(p, q, repeats):
#Setup an empty array
KL_scores = np.zeros(repeats)
#Now repeat the KL calulation of the adjusted arrays
for i in range(repeats):
#Match the array size
pMatched = uneven_arrays(p,q)[0]
qMatched = uneven_arrays(p,q)[1]
#Calculate KL divergence of matched arrays, in base 2 since it is nicer to work with
KL_scores[i] = entropy(pMatched, qMatched, base=2)
#Remove any unusually large values from the score
KL_scores = KL_scores[KL_scores < 1e2]
#Return the average KL divergence amongst these
return np.mean(KL_scores)
Note: KL(p||q) asks how different an estimate distribution q is from the reference distribution p. As AMIP is our reference, we take p to this, and q to be the historical output.
KL_score = uneven_KL(p=amip_gaussian, q=historical_gaussian,repeats=10000)
print('The GISS model achieved a KL score of:', np.round(KL_score,3))
The GISS model achieved a KL score of: 0.853
def KL_closed_form(mean_1, cov_1, mean_2, cov_2):
mean_1 = np.reshape(mean_1, (3, 1))
mean_2 = np.reshape(mean_2, (3, 1))
first_part = np.log(
np.linalg.det(cov_2) / np.linalg.det(cov_1)) + np.trace(
np.dot((np.linalg.inv(cov_2)), cov_1))
second_part = np.dot(
np.dot((np.transpose(mean_2 - mean_1)), np.linalg.inv(cov_2)),
(mean_2 - mean_1)) - len(mean_1)
KL = 0.5 * (first_part + second_part)
return KL[0,0]
KL_closed_form(amip_mean, amip_cov, historical_mean, historical_cov)
6.664769474125807
def uneven_JS(p,q,repeats):
#Setup an empty array
JS_scores = np.zeros(repeats)
#Now repeat the JS calulation of the adjusted arrays
for i in range(repeats):
#Match the array size
pMatched = uneven_arrays(p,q)[0]
qMatched = uneven_arrays(p,q)[1]
#Calculate JS of matched arrays
JS_scores[i] = jensenshannon(pMatched, qMatched)
#Remove any unusually large values from the score
JS_scores = JS_scores[JS_scores < 1e2]
#Return the average KL divergence amongst these
return np.mean(JS_scores)
JS_score = uneven_JS(p=amip_gaussian,
q=historical_gaussian,
repeats=10000)
print('The GISS model achieved a JS score of:', np.round(JS_score,3))
The GISS model achieved a JS score of: 0.346
To calculate the Wasserstein distance, we ignore the fitted Gaussian and perform the following steps:
(1): Set the 3D points in each experiment to have a uniformly distributed mass
(2): Calculate the cost matrix associated to the position of the 3D points in the two experiments
(3): Use the weights and cost matrix to calculate the (approximate) Wasserstein distance/optimal transport cost.
amip_weights = np.ones((len(amip_outputs[:,0]), )) / len(amip_outputs[:,0])
historical_weights = np.ones((len(historical_outputs[:,0]), )) / len(historical_outputs[:,0])
%%time
#Cheap fix to stop kernel dying:
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
#Generate the cost matrix with p=1
cost_matrix = ot.dist(amip_outputs, historical_outputs, metric='minkowski', p=1)
#Divide by the maximum value
cost_matrix /= cost_matrix.max()
#Plot the cost matrix:
plt.imshow(cost_matrix, interpolation ='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.5)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('Cost matrix')
#plt.savefig('images/TasPrPsl_cost_matrix_p=1', bbox_inches='tight', dpi=300)
#Now compute Wasserstein distance, which takes the weights and the cost matrix
#NOTE: WE MUST TAKE THE pTH ROOT IN ORDER TO GET THE REAL WASSERSTEIN DISTANCE
wasserstein1_score = ot.emd2(amip_weights, historical_weights,
cost_matrix)**(1/1)
print('The GISS model achieved a Wasserstein-1 score of:', np.round(wasserstein1_score,3))
The GISS model achieved a Wasserstein-1 score of: 0.371 Wall time: 37.9 ms
#Compute the optimal transport matrix:
optimal_transport = ot.emd(amip_weights, historical_weights,
cost_matrix)
plt.imshow(optimal_transport, interpolation='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.7)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('GISS Optimal Transport Matrix')
plt.savefig('images/TasPrPsl_optimal_transport_matrix_p=1', bbox_inches='tight', dpi=300)
plt.figure(1, figsize=(15, 4))
plt.suptitle("GISS Optimal Transport Paths", fontsize="xx-large")
#Tas vs Pr
plt.subplot(1, 3, 1)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,1]], historical_outputs[:,[0,1]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 1], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 1], 'ob', label='Historical samples')
plt.title('$p=1$', rotation='vertical', x=-0.25,y=0.4, fontsize='xx-large')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
#Tas vs Psl
plt.subplot(1, 3, 2)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,2]], historical_outputs[:,[0,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#Pr vs Psl
plt.subplot(1, 3, 3)
ot.plot.plot2D_samples_mat(amip_outputs[:,[1,2]], historical_outputs[:,[1,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 1], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 1], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#plt.legend(loc='upper left')
plt.tight_layout()
plt.savefig('images/TasPrPsl_optimal_transport_paths_p=1.eps', format='eps', bbox_inches='tight', dpi=300)
The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.
%%time
#Cheap fix to stop kernel dying:
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
#Generate the cost matrix with p=2
cost_matrix = ot.dist(amip_outputs, historical_outputs, metric='minkowski', p=2)
#Divide by the maximum value
cost_matrix /= cost_matrix.max()
#Plot the cost matrix:
plt.imshow(cost_matrix, interpolation ='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.5)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('Cost matrix')
plt.savefig('images/TasPrPsl_cost_matrix_p=2', bbox_inches='tight', dpi=300)
#Now compute Wasserstein distance, which takes the weights and the cost matrix
#NOTE: WE MUST TAKE THE pTH ROOT IN ORDER TO GET THE REAL WASSERSTEIN DISTANCE
wasserstein2_score = ot.emd2(amip_weights, historical_weights,
cost_matrix)**(1/2)
print('The GISS model achieved a Wasserstein-2 score of:', np.round(wasserstein2_score,3))
The GISS model achieved a Wasserstein-2 score of: 0.573 Wall time: 308 ms
#Compute the optimal transport matrix:
optimal_transport = ot.emd(amip_weights, historical_weights,
cost_matrix)
plt.imshow(optimal_transport, interpolation='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.7)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('GISS Optimal Transport Matrix')
plt.savefig('images/TasPrPsl_optimal_transport_matrix_p=2', bbox_inches='tight', dpi=300)
plt.figure(1, figsize=(15, 4))
#Tas vs Pr
plt.subplot(1, 3, 1)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,1]], historical_outputs[:,[0,1]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 1], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 1], 'ob', label='Historical samples')
plt.title('$p=2$', rotation='vertical', x=-0.25,y=0.4, fontsize='xx-large')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
#Tas vs Psl
plt.subplot(1, 3, 2)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,2]], historical_outputs[:,[0,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#Pr vs Psl
plt.subplot(1, 3, 3)
ot.plot.plot2D_samples_mat(amip_outputs[:,[1,2]], historical_outputs[:,[1,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 1], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 1], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.tight_layout()
plt.savefig('images/TasPrPsl_optimal_transport_paths_p=2.eps', format='eps', bbox_inches='tight', dpi=300)
The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.
%%time
#Cheap fix to stop kernel dying:
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
#Generate the cost matrix with p=3
cost_matrix = ot.dist(amip_outputs, historical_outputs, metric='minkowski', p=3)
#Divide by the maximum value
cost_matrix /= cost_matrix.max()
#Plot the cost matrix:
plt.imshow(cost_matrix, interpolation ='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.5)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('Cost matrix')
plt.savefig('images/TasPrPsl_cost_matrix_p=3', bbox_inches='tight', dpi=300)
#Now compute Wasserstein distance, which takes the weights and the cost matrix
#NOTE: WE MUST TAKE THE pTH ROOT IN ORDER TO GET THE REAL WASSERSTEIN DISTANCE
wasserstein3_score = ot.emd2(amip_weights, historical_weights,
cost_matrix)**(1/3)
print('The GISS model achieved a Wasserstein-3 score of:', np.round(wasserstein3_score,3))
The GISS model achieved a Wasserstein-3 score of: 0.685 Wall time: 306 ms
#Compute the optimal transport matrix:
optimal_transport = ot.emd(amip_weights, historical_weights,
cost_matrix)
plt.imshow(optimal_transport, interpolation='nearest')
#Add a colorbar below
cbar = plt.colorbar(shrink=0.7)
cbar.set_label('Cost', rotation=270,labelpad=15)
plt.xlabel('Historical Ensemble Member')
plt.ylabel('AMIP Ensemble Member')
plt.title('GISS Optimal Transport Matrix')
plt.savefig('images/TasPrPsl_optimal_transport_matrix_p=3', bbox_inches='tight', dpi=300)
plt.figure(1, figsize=(15, 4))
#Tas vs Pr
plt.subplot(1, 3, 1)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,1]], historical_outputs[:,[0,1]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 1], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 1], 'ob', label='Historical samples')
plt.title('$p=3$', rotation='vertical', x=-0.25,y=0.4, fontsize='xx-large')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
#Tas vs Psl
plt.subplot(1, 3, 2)
ot.plot.plot2D_samples_mat(amip_outputs[:,[0,2]], historical_outputs[:,[0,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 0], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 0], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Tas (K)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#Pr vs Psl
plt.subplot(1, 3, 3)
ot.plot.plot2D_samples_mat(amip_outputs[:,[1,2]], historical_outputs[:,[1,2]], optimal_transport, c=[.7, .7, 1])
plt.plot(amip_outputs[:, 1], amip_outputs[:, 2], 'r^', label='AMIP samples')
plt.plot(historical_outputs[:, 1], historical_outputs[:, 2], 'ob', label='Historical samples')
plt.xlabel('Pr ($g m^{-2} d^{-1}$)', fontsize='x-large')
plt.ylabel('Psl (Pa)', fontsize='x-large')
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.tight_layout()
plt.savefig('images/TasPrPsl_optimal_transport_paths_p=3.eps', format='eps', bbox_inches='tight', dpi=300)
The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque. The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.
wasserstein_scores = [wasserstein1_score, wasserstein2_score, wasserstein3_score]
np.round(wasserstein_scores,3)
array([0.371, 0.573, 0.685])
My implementation, using the empirical estimators of expectation found by summing pairwise differences:
def calculate_energy_statistics(x, y):
n = len(x)
m = len(y)
A = 0
B = 0
C = 0
for i in range(n):
for j in range(m):
A += np.linalg.norm(x[i] - y[j])
A = A / (n * m)
for i in range(n):
for j in range(n):
B += np.linalg.norm(x[i] - x[j])
B = B / (n * n)
for i in range(m):
for j in range(m):
C += np.linalg.norm(y[i] - y[j])
C = C / (m * m)
energy_distance = np.sqrt(2 * A - B - C)
test_statistic = (2 * A - B - C) / (2 * A)
return {
'Energy distance': np.round(energy_distance, 3),
'Test statistic': np.round(test_statistic, 4)
}
calculate_energy_statistics(amip_outputs, historical_outputs)
{'Energy distance': 14.987, 'Test statistic': 0.4901}